setwd("/Users/yaping/Documents/lab/project/bisulfite-snp/roc/roc20110804")
library(ROC)
##make methy cytosine caller roc
bs_cytosine<-read.table("vcf_methy_summary.all.anotherway.txt",header=F,sep="\t")
bismark_cytosine<-read.table("bismark_methy_summary.txt",header=F,sep="\t")
berman_paper<-read.table("normal_berman_paper_methy_summary.txt",header=F,sep="\t")
#berman_withoutref<-read.table("normal_berman_without_refc_methy_summary.txt",header=F,sep="\t")
temp<-bismark_cytosine[order(bismark_cytosine[,3]),]
x<-1-temp[,3]
y<-1-temp[,4]
rocBsAllBismark<-new('rocc',sens=y,spec=x)

plot(rocBsAllBismark,col="blue",xlim=c(0,0.4),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,axes=F,xlab="",ylab="")

par(new=T)

temp<-berman_paper[order(berman_paper[,7]),]
x<-1-temp[,7]
y<-1-temp[,8]
rocBerman_paper<-new('rocc',sens=y,spec=x)

plot(rocBerman_paper,col="green",xlim=c(0,0.4),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,axes=F,xlab="",ylab="")

#par(new=T)

#temp<-berman_withoutref[order(berman_withoutref[,7]),]
#x<-1-temp[,7]
#y<-1-temp[,8]
#rocBerman_withoutref<-new('rocc',sens=y,spec=x)

#plot(rocBerman_withoutref,col="yellow",xlim=c(0,0.4),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,axes=F,xlab="",ylab="")

par(new=T)

temp<-bs_cytosine[order(bs_cytosine[,8]),]
x<-1-temp[,8]
y<-1-temp[,11]
rocBsAll<-new('rocc',sens=y,spec=x)

plot(rocBsAll,col="red",xlim=c(0,0.4),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,xlab="False Positive Rate",ylab="1-False Negative Rate", main="ROC curve on polymorphic cytosines (1M SNP array) ")


legend("bottomright",c("Bis-SNP","Bismark","2011_Berman"),col=c("red","blue","green"),lty=1,lwd=2,cex=1.2)

point_bismark<-bismark_cytosine[bismark_cytosine[,2]==0,]
points(point_bismark[,3], 1-point_bismark[,4],type="p",col="blue",pch=21,cex=1.5)
text(point_bismark[,3]+0.013, 1-point_bismark[,4],"y",col="blue",font=2,cex=1.5)

point_berman<-berman_paper[berman_paper[,5]<0.91 & berman_paper[,5]>0.895,]
points(point_berman[,7], 1-point_berman[,8],type="p",col="green",pch=21,cex=1.5)
text(point_berman[,7]+0.013, 1-point_berman[,8]-0.02,"x",col="green",font=2,cex=1.5)

#point_bissnp<-bs_cytosine[bs_cytosine[,5]==20,]
#points(point_bissnp[,8], 1-point_bissnp[,11],type="p",col="red",pch=21,cex=1.5)
#text(point_bissnp[,8]+0.013, 1-point_bissnp[,11],"a",col="red",font=2,cex=1.5)

tempCut1<-berman_paper[berman_paper[,5]<0.91 & berman_paper[,5]>0.895,]
abline(v=tempCut1[,7])
text(tempCut1[,7],0.2,"minCTcontent>0.90 & minGoppositeStrand>0.90")

tempCut2<-bs_cytosine[bs_cytosine[,5]==20,]
abline(v=tempCut2[,8],col="purple")
text(tempCut2[,8],0.4,"Bis-SNP_confidance>=20",col="purple")

tempCut3<-bs_cytosine[bs_cytosine[,5]==10,]
abline(v=tempCut3[,8],col="purple")
text(tempCut3[,8],0.6,"Bis-SNP_confidance>=10",col="purple")

#zoom in

temp<-bismark_cytosine[order(bismark_cytosine[,3]),]
x<-1-temp[,3]
y<-1-temp[,4]
rocBsAllBismark<-new('rocc',sens=y,spec=x)

plot(rocBsAllBismark,col="blue",xlim=c(0,0.027),ylim=c(0.7,0.8),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,axes=F,xlab="",ylab="")

par(new=T)

temp<-berman_paper[order(berman_paper[,7]),]
x<-1-temp[,7]
y<-1-temp[,8]
rocBerman_paper<-new('rocc',sens=y,spec=x)

plot(rocBerman_paper,col="green",xlim=c(0,0.027),ylim=c(0.7,0.8),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,axes=F,xlab="",ylab="")

#par(new=T)

#temp<-berman_withoutref[order(berman_withoutref[,7]),]
#x<-1-temp[,7]
#y<-1-temp[,8]
#rocBerman_withoutref<-new('rocc',sens=y,spec=x)

#plot(rocBerman_withoutref,col="yellow",xlim=c(0,0.027),ylim=c(0.7,0.8),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,axes=F,xlab="",ylab="")

par(new=T)

temp<-bs_cytosine[order(bs_cytosine[,8]),]
x<-1-temp[,8]
y<-1-temp[,11]
rocBsAll<-new('rocc',sens=y,spec=x)

plot(rocBsAll,col="red",xlim=c(0,0.027),ylim=c(0.7,0.8),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,xlab="False Positive Rate",ylab="1-False Negative Rate", main="ROC curve on polymorphic cytosines (zoom in) ")


legend("bottomright",c("Bis-SNP","Bismark","2011_Berman"),col=c("red","blue","green"),lty=1,lwd=2,cex=1.2)

tempCut1<-berman_paper[berman_paper[,5]<0.91 & berman_paper[,5]>0.895,]
abline(v=tempCut1[,7])
text(tempCut1[,7]+0.002,0.73,"minCTcontent>0.90 & minGoppositeStrand>0.90")

tempCut2<-bs_cytosine[bs_cytosine[,5]==20,]
abline(v=tempCut2[,8],col="purple")
text(tempCut2[,8],0.76,"Bis-SNP_confidance>=20",col="purple")

tempCut3<-bs_cytosine[bs_cytosine[,5]==10,]
abline(v=tempCut3[,8],col="purple")
text(tempCut3[,8],0.78,"Bis-SNP_confidance>=10",col="purple")

##make snp roc

bs_all<-read.table("vcf_snp_summary.all.txt",header=F,sep="\t")
kallele_all<-read.table("normal_kallele_snp_summary.txt",header=F,sep="\t")
pash_all<-read.table("vcf_PASH_snp_summary.txt",header=F,sep="\t")
temp<-kallele_all[order(kallele_all[,8]),]
x<-1-temp[,8]
y<-1-temp[,9]
rocKalleleAll<-new('rocc',sens=y,spec=x)
plot(rocKalleleAll,col="green",xlim=c(0,0.0065),ylim=c(0,1),axes=F,xlab="",ylab="",lty=1,font=2,lwd=2)



par(new=T)
temp<-bs_all[order(bs_all[,8]),]
x<-1-temp[,8]
y<-1-temp[,11]
rocBsAll<-new('rocc',sens=y,spec=x)

plot(rocBsAll,col="red",xlim=c(0,0.0065),ylim=c(0,1),axes=F,xlab="",ylab="",lty=1,font=2,lwd=2)

par(new=T)
temp<-pash_all[order(pash_all[,4]),]
x<-1-temp[,4]
y<-1-temp[,5]
rocPashAll<-new('rocc',sens=y,spec=x)

plot(rocPashAll,col="blue",xlim=c(0,0.0065),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,xlab="False Positive Rate",ylab="1-False Negative Rate", main="ROC curve on heterozygous SNP sites (1M SNP array)")

legend("bottomright",c("Bis-SNP","K-allele","2010_Shoemaker & Samtools"),col=c("red","green","blue"),lty=1,lwd=2,cex=1.2)

tempCut1<-kallele_all[kallele_all[,5]==2,]
abline(v=tempCut1[,8])
text(tempCut1[,8]+0.0016,0.2,"AlleleCount>=2 & AlleleFreq>=0.2")
tempCut2<-kallele_all[kallele_all[,5]==1,]
#abline(v=tempCut2[,8])
#text(tempCut2[,8]+0.005,0.2,"AlleleCount>=1 & AlleleFreq>=0.1")

tempCut3<-kallele_all[kallele_all[,5]==3,]
abline(v=tempCut3[,8])
text(tempCut3[,8]+0.0016,0.4,"AlleleCount>=3 & AlleleFreq>=0.3")

tempCut4<-bs_all[bs_all[,5]==20,]
abline(v=tempCut4[,8],col="purple")
text(tempCut4[,8],0.6,"Bis-SNP_confidance>=20",col="purple")

tempCut5<-bs_all[bs_all[,5]==10,]
abline(v=tempCut5[,8],col="purple")
text(tempCut5[,8],0.8,"Bis-SNP_confidance>=10",col="purple")


#zoom in 
temp<-kallele_all[order(kallele_all[,8]),]
x<-1-temp[,8]
y<-1-temp[,9]
rocKalleleAll<-new('rocc',sens=y,spec=x)

plot(rocKalleleAll,col="green",xlim=c(0,tempCut3[,8]),ylim=c(0.0,0.2),axes=F,xlab="",ylab="",lty=1,font=2,lwd=2)

par(new=T)
temp<-pash_all[order(pash_all[,4]),]
x<-1-temp[,4]
y<-1-temp[,5]
rocPashAll<-new('rocc',sens=y,spec=x)

plot(rocPashAll,col="blue",xlim=c(0,tempCut3[,8]),ylim=c(0.0,0.2),axes=F,xlab="",ylab="",lty=1,font=2,lwd=2)

par(new=T)
temp<-bs_all[order(bs_all[,8]),]
x<-1-temp[,8]
y<-1-temp[,11]
rocBsAll<-new('rocc',sens=y,spec=x)

plot(rocBsAll,col="red",xlim=c(0,tempCut3[,8]),ylim=c(0.0,0.2),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,xlab="False Positive Rate",ylab="1-False Negative Rate", main="ROC curve on heterozygous SNP sites (zoom in)")
legend("bottomright",c("Bis-SNP","K-allele","2010_Coarfa & samtools"),col=c("red","green","blue"),lty=1,lwd=2,cex=1.2)
tempCut1<-kallele_all[kallele_all[,5]==2,]
abline(v=tempCut1[,8])
text(tempCut1[,8]-0.00053,0.8,"AlleleCount>=2 & AlleleFreq>=0.2")
#tempCut2<-kallele_all[kallele_all[,5]==1,]
#abline(v=tempCut2[,8])
#text(tempCut2[,8]-0.0035,0.90,"AlleleCount>=1 & AlleleFreq>=0.1")

tempCut3<-kallele_all[kallele_all[,5]==3,]
abline(v=tempCut3[,8])
text(tempCut3[,8]-0.00053,0.10,"AlleleCount>=3 & AlleleFreq>=0.3")



##make ctsnp
bs_ctsnp<-read.table("vcf_snp_summary.ctsnp.txt",header=F,sep="\t")
kallele_ctsnp<-read.table("normal_kallele_snp_summary.ctsnp.txt",header=F,sep="\t")
#pash_all<-read.table("vcf_PASH_snp_summary.txt",header=F,sep="\t")
temp<-kallele_ctsnp[order(kallele_ctsnp[,8]),]
x<-1-temp[,8]
y<-1-temp[,9]
rocKallelectsnp<-new('rocc',sens=y,spec=x)
plot(rocKallelectsnp,col="green",xlim=c(0,0.0065),ylim=c(0,1),axes=F,xlab="",ylab="",lty=1,font=2,lwd=2)

#par(new=T)
#temp<-pash_all[order(pash_all[,4]),]
#x<-1-temp[,4]
#y<-1-temp[,5]
#rocPashAll<-new('rocc',sens=y,spec=x)

#plot(rocPashAll,col="blue",xlim=c(0,0.03),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,axes=F,xlab="",ylab="")

par(new=T)
temp<-bs_ctsnp[order(bs_ctsnp[,8]),]
x<-1-temp[,8]
y<-1-temp[,11]
rocBsctsnp<-new('rocc',sens=y,spec=x)

plot(rocBsctsnp,col="red",xlim=c(0,0.0065),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,xlab="False Positive Rate",ylab="1-False Negative Rate", main="ROC curve on heterozygous CT SNP sites")

legend("bottomright",c("Bis-SNP","K-allele"),col=c("red","green"),lty=1,lwd=2,cex=1.2)

tempCut1<-kallele_all[kallele_all[,5]==2,]
abline(v=tempCut1[,8])
text(tempCut1[,8]+0.0016,0.2,"AlleleCount>=2 & AlleleFreq>=0.2")
tempCut2<-kallele_all[kallele_all[,5]==1,]
#abline(v=tempCut2[,8])
#text(tempCut2[,8]+0.005,0.2,"AlleleCount>=1 & AlleleFreq>=0.1")

tempCut3<-kallele_all[kallele_all[,5]==3,]
abline(v=tempCut3[,8])
text(tempCut3[,8]+0.0016,0.4,"AlleleCount>=3 & AlleleFreq>=0.3")


#zoom in 
temp<-kallele_ctsnp[order(kallele_ctsnp[,8]),]
x<-1-temp[,8]
y<-1-temp[,9]
rocKallelectsnp<-new('rocc',sens=y,spec=x)

plot(rocKallelectsnp,col="green",xlim=c(0,tempCut1[,8]),ylim=c(0.6,1),axes=F,xlab="",ylab="",lty=1,font=2,lwd=2)

par(new=T)
temp<-bs_ctsnp[order(bs_ctsnp[,8]),]
x<-1-temp[,8]
y<-1-temp[,11]
rocBsctsnp<-new('rocc',sens=y,spec=x)

plot(rocBsctsnp,col="red",xlim=c(0,tempCut1[,8]),ylim=c(0.6,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,xlab="False Positive Rate",ylab="1-False Negative Rate", main="ROC curve on heterozygous CT SNP sites (zoom in)")
legend("bottomright",c("Bis-SNP","K-allele"),col=c("red","green"),lty=1,lwd=2,cex=1.2)
tempCut1<-kallele_all[kallele_all[,5]==2,]
abline(v=tempCut1[,8])
text(tempCut1[,8]-0.00053,0.8,"AlleleCount>=2 & AlleleFreq>=0.2")
#tempCut2<-kallele_all[kallele_all[,5]==1,]
#abline(v=tempCut2[,8])
#text(tempCut2[,8]-0.0035,0.90,"AlleleCount>=1 & AlleleFreq>=0.1")

tempCut3<-kallele_all[kallele_all[,5]==3,]
abline(v=tempCut3[,8])
text(tempCut3[,8]+0.00053,0.70,"AlleleCount>=3 & AlleleFreq>=0.3")



##make nonctsnp
bs_nonctsnp<-read.table("vcf_snp_summary.nonctsnp.txt",header=F,sep="\t")
kallele_nonctsnp<-read.table("normal_kallele_snp_summary.nonctsnp.txt",header=F,sep="\t")
#pash_all<-read.table("vcf_PASH_snp_summary.txt",header=F,sep="\t")
temp<-kallele_nonctsnp[order(kallele_nonctsnp[,8]),]
x<-1-temp[,8]
y<-1-temp[,9]
rocKallelenonctsnp<-new('rocc',sens=y,spec=x)
plot(rocKallelenonctsnp,col="green",xlim=c(0,0.0065),ylim=c(0,1),axes=F,xlab="",ylab="",lty=1,font=2,lwd=2)

#par(new=T)
#temp<-pash_all[order(pash_all[,4]),]
#x<-1-temp[,4]
#y<-1-temp[,5]
#rocPashAll<-new('rocc',sens=y,spec=x)

#plot(rocPashAll,col="blue",xlim=c(0,0.0065),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,axes=F,xlab="",ylab="")

par(new=T)
temp<-bs_nonctsnp[order(bs_nonctsnp[,8]),]
x<-1-temp[,8]
y<-1-temp[,11]
rocBsnonctsnp<-new('rocc',sens=y,spec=x)

plot(rocBsnonctsnp,col="red",xlim=c(0,0.0065),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,xlab="False Positive Rate",ylab="1-False Negative Rate", main="ROC curve on heterozygous Non-CT SNP sites")

legend("bottomright",c("Bis-SNP","K-allele"),col=c("red","green"),lty=1,lwd=2,cex=1.2)

tempCut1<-kallele_all[kallele_all[,5]==2,]
abline(v=tempCut1[,8])
text(tempCut1[,8]+0.0016,0.2,"AlleleCount>=2 & AlleleFreq>=0.2")
tempCut2<-kallele_all[kallele_all[,5]==1,]
#abline(v=tempCut2[,8])
#text(tempCut2[,8]+0.005,0.2,"AlleleCount>=1 & AlleleFreq>=0.1")

tempCut3<-kallele_all[kallele_all[,5]==3,]
abline(v=tempCut3[,8])
text(tempCut3[,8]+0.0016,0.4,"AlleleCount>=3 & AlleleFreq>=0.3")


#zoom in 
temp<-kallele_nonctsnp[order(kallele_nonctsnp[,8]),]
x<-1-temp[,8]
y<-1-temp[,9]
rocKallelenonctsnp<-new('rocc',sens=y,spec=x)

plot(rocKallelenonctsnp,col="green",xlim=c(0,tempCut1[,8]),ylim=c(0.6,1),axes=F,xlab="",ylab="",lty=1,font=2,lwd=2)

par(new=T)
temp<-bs_nonctsnp[order(bs_nonctsnp[,8]),]
x<-1-temp[,8]
y<-1-temp[,11]
rocBsnonctsnp<-new('rocc',sens=y,spec=x)

plot(rocBsnonctsnp,col="red",xlim=c(0,tempCut1[,8]),ylim=c(0.6,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,xlab="False Positive Rate",ylab="1-False Negative Rate", main="ROC curve on heterozygous Non-CT SNP sites (zoom in)")
legend("bottomright",c("Bis-SNP","K-allele"),col=c("red","green"),lty=1,lwd=2,cex=1.2)
tempCut1<-kallele_all[kallele_all[,5]==2,]
abline(v=tempCut1[,8])
text(tempCut1[,8]-0.00053,0.8,"AlleleCount>=2 & AlleleFreq>=0.2")
#tempCut2<-kallele_all[kallele_all[,5]==1,]
#abline(v=tempCut2[,8])
#text(tempCut2[,8]-0.0035,0.90,"AlleleCount>=1 & AlleleFreq>=0.1")

tempCut3<-kallele_all[kallele_all[,5]==3,]
abline(v=tempCut3[,8])
text(tempCut3[,8]+0.00053,0.70,"AlleleCount>=3 & AlleleFreq>=0.3")




###################
#use ROCR lib
lib(ROCR)


bs_all<-read.table("vcf_snp_summary.all.txt",header=F,sep="\t")
kallele_all<-read.table("normal_kallele_snp_summary.txt",header=F,sep="\t")
pash_all<-read.table("vcf_PASH_snp_summary.txt",header=F,sep="\t")
temp<-kallele_all[order(kallele_all[,8]),]
x<-temp[,8]
y<-1-temp[,9]
rocKalleleAll<-new("performance",x.values= list(x),y.values= list(y), alpha.values = list(temp[,5]),x.name="False Positive",y.name="True Positive",alpha.name = "Confidance level")
plot(rocKalleleAll,col="green",xlim=c(0,0.03),ylim=c(0,1),axes=F,xlab="",ylab="",lty=1,font=2,lwd=2)

#par(new=T)
#temp<-pash_all[order(pash_all[,4]),]
#x<-1-temp[,4]
#y<-1-temp[,5]
#rocPashAll<-new('rocc',sens=y,spec=x)

#plot(rocPashAll,col="blue",xlim=c(0,0.03),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,axes=F,xlab="",ylab="")

par(new=T)
temp<-bs_all[order(bs_all[,8]),]
x<-temp[,8]
y<-1-temp[,11]
rocBsAll<-new("performance",x.values= list(x),y.values= list(y), alpha.values = list(temp[,5]),x.name="False Positive",y.name="True Positive",alpha.name = "Confidance level")

plot(rocBsAll,col="red",xlim=c(0,0.03),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,xlab="False Positive Rate",ylab="1-False Negative Rate", main="ROC curve on heterozygous SNP sites")

legend("bottomright",c("Bis-SNP","K-allele"),col=c("red","green"),lty=1,lwd=2,cex=1.2)

tempCut1<-kallele_all[kallele_all[,5]==2,]
abline(v=tempCut1[,8])
text(tempCut1[,8]+0.007,0.4,"AlleleCount>=2 & AlleleFreq>=0.2")
tempCut2<-kallele_all[kallele_all[,5]==1,]
abline(v=tempCut2[,8])
text(tempCut2[,8]+0.005,0.2,"AlleleCount>=1 & AlleleFreq>=0.1")

tempCut3<-kallele_all[kallele_all[,5]==3,]
abline(v=tempCut3[,8])
text(tempCut3[,8]+0.007,0.6,"AlleleCount>=3 & AlleleFreq>=0.3")





###make ROC curve in different downsampling
setwd("/Users/yaping/Documents/lab/project/bisulfite-snp/roc/roc20110826")
library(ROC)
bs_all<-read.table("vcf_snp_summary.all.downsampling.txt",header=F,sep="\t")
legName<-c()
legPch<-c()
for(i in c(1,5,10,15,20,25,30,35,40)){
	temp<-bs_all[bs_all[,16]==i,]
	temp<-temp[order(temp[,8]),]
	x<-1-temp[,8]
	y<-1-temp[,11]
	rocBsAll<-new('rocc',sens=y,spec=x)
	plot(temp[,8],y,xlim=c(0,0.01),ylim=c(0,1),axes=F,xlab="",ylab="",lty=1,font=2,lwd=2,pch=i/5,type="b")
#	plot(rocBsAll,xlim=c(0,0.01),ylim=c(0,1),axes=F,xlab="",ylab="",lty=1,font=2,lwd=2,pch=i)
	name<-paste("read_cov <= ",i,sep="")
	legName<-c(legName,name)
	legPch<-c(legPch,i/5)
	par(new=T)
}


#plot(rocPashAll,col="blue",xlim=c(0,1),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,xlab="False Positive Rate",ylab="1-False Negative Rate", main="ROC curve on heterozygous SNP sites (1M SNP array)")
axis(1,at=seq(0,0.01,by=0.002),pch=16,lty=1,font=2,cex.axis=1.2)
axis(2,at=seq(0,1,by=0.1),pch=16,lty=1,font=2,cex.axis=1.2)
legend("bottomright",legName,pch=legPch,lty=1,lwd=2,cex=1.2)




##make for jone lab meeting
setwd("/Users/yaping/Documents/lab/project/bisulfite-snp/roc/roc20110810")
bs_cytosine<-read.table("vcf_methy_summary.all.anotherway.txt",header=F,sep="\t")
bismark_cytosine<-read.table("bismark_methy_summary.txt",header=F,sep="\t")
berman_paper<-read.table("normal_berman_paper_methy_summary.txt",header=F,sep="\t")
#berman_withoutref<-read.table("normal_berman_without_refc_methy_summary.txt",header=F,sep="\t")
temp<-bismark_cytosine[order(bismark_cytosine[,3]),]
x<-1-temp[,3]
y<-1-temp[,4]
rocBsAllBismark<-new('rocc',sens=y,spec=x)

plot(rocBsAllBismark,col="blue",xlim=c(0,0.1),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,axes=F,xlab="",ylab="")

par(new=T)

temp<-berman_paper[order(berman_paper[,7]),]
x<-1-temp[,7]
x<-c(1,x)
y<-1-temp[,8]
y<-c(0,y)
rocBerman_paper<-new('rocc',sens=y,spec=x)

plot(rocBerman_paper,col="purple",xlim=c(0,0.1),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,axes=F,xlab="",ylab="")

#par(new=T)

#temp<-berman_withoutref[order(berman_withoutref[,7]),]
#x<-1-temp[,7]
#y<-1-temp[,8]
#rocBerman_withoutref<-new('rocc',sens=y,spec=x)

#plot(rocBerman_withoutref,col="yellow",xlim=c(0,0.1),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,axes=F,xlab="",ylab="")

par(new=T)

temp<-bs_cytosine[order(bs_cytosine[,8]),]
x<-1-temp[,8]
y<-1-temp[,11]
rocBsAll<-new('rocc',sens=y,spec=x)

plot(rocBsAll,col="red",xlim=c(0,0.1),ylim=c(0,1),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,xlab="False Positive Rate",ylab="1-False Negative Rate", main="ROC curve on polymorphic cytosines (1M SNP array) ")


legend("bottomright",c("Bis-SNP","Bismark","2011_Berman"),col=c("red","blue","purple"),lty=1,lwd=2,cex=1.2)
